Data Preparation

# smacof analysis of 1988 wine data
# Lew Harvey
# 3 February 2016

library("smacof")
library("rgl")

Read in the data and wine properties

fn <- "wine_88.csv"
df <- read.csv(fn, header = TRUE)
df$w10 <- as.numeric(df$w10)

# read in wine properties
df.prop <- read.csv("wine_properties.csv", header = TRUE)
df.prop$wine <- as.character(df.prop$wine)
df.prop$color <- as.character(df.prop$color)

The data in the input file looks like this for the first subject:

head(df, 10)
##    w01 w02 w03 w04 w05 w06 w07 w08 w09 w10
## 1    0  NA  NA  NA  NA  NA  NA  NA  NA  NA
## 2    2   0  NA  NA  NA  NA  NA  NA  NA  NA
## 3    2   2   0  NA  NA  NA  NA  NA  NA  NA
## 4    4   4   2   0  NA  NA  NA  NA  NA  NA
## 5    2   2   6   2   0  NA  NA  NA  NA  NA
## 6    5   4   6   4   7   0  NA  NA  NA  NA
## 7    3   5   7   7   5   3   0  NA  NA  NA
## 8    5   6   7   6   5   3   7   0  NA  NA
## 9    3   5   6   6   4   2   6   5   0  NA
## 10   2   4   3   3   7   5   6   6   2   0

Convert the data into a list of matrices, one for each subject. This is a bit more complicated than it should be, but the original data are in an SPSS file that requires some massaging to get the data into a form for R

# loop through the data frame to
# create separate matrices for the
# data sets of each subject

nstim <- ncol(df)
nsets <- nrow(df) / ncol(df)
mlist <- list(NA)
mm <- list(NA)
winename <- colnames(df)

for (i in 1:nsets){
    mm[i] <- paste("m", i, sep = "")
    i1 <- ((i - 1) * nstim) + 1
    i2 <- i1 + nstim - 1
    mlist[i] <- list(mm = as.matrix(df[i1:i2,]))
    row.names(mlist[[i]]) <- winename
    mlist[[i]] <- as.dist(mlist[[i]])
}
names(mlist) <- unlist(mm)

Now the matrices are in a list, one for each subject. Here are the first subject’s data as a distance matrix:

mlist[1]
## $m1
##     w01 w02 w03 w04 w05 w06 w07 w08 w09
## w02   2                                
## w03   2   2                            
## w04   4   4   2                        
## w05   2   2   6   2                    
## w06   5   4   6   4   7                
## w07   3   5   7   7   5   3            
## w08   5   6   7   6   5   3   7        
## w09   3   5   6   6   4   2   6   5    
## w10   2   4   3   3   7   5   6   6   2

2D INDSCAL Solution

wine2D <- smacofIndDiff(mlist,
                      ndim = 2,
                      type = "ordinal",
                      constraint = "indscal",
                      itmax = 1000)
wine2D
## 
## Call: smacofIndDiff(delta = mlist, ndim = 2, type = "ordinal", constraint = "indscal", 
##     itmax = 1000)
## 
## Model: Three-way SMACOF 
## Number of objects: 10 
## Stress-1 value: 0.2201172 
## Number of iterations: 440

Plot 2D Solution and Fits

par(mfcol = c(2, 2))
par(pty = "s")
plot(wine2D, type = "p",
     plot.type = "confplot",
     label.conf = c(TRUE, 1, 1),
     xlim = c(-0.9, 0.9),
     ylim = c(-0.9, 0.9))
abline(h = 0, col = "gray")
abline(v = 0, col = "gray")
par(pty = "m")

plot(wine2D, plot.type = "Shepard") 
plot(wine2D, plot.type = "stressplot", ylim = c(0, 25))
plot(wine2D, plot.type = "resplot")

par(mfcol = c(1, 1))

Configuration Plot without Color

par(pty = "s")
plot(wine2D, type = "p", pch = 19,
     plot.type = "confplot",
     label.conf = c(TRUE, 1, 1),
     xlim = c(-0.9, 0.9),
     ylim = c(-0.9, 0.9))
abline(h = 0, col = "gray")
abline(v = 0, col = "gray")

par(pty = "m")

Configuration Plot with Color

par(pty = "s")
plot(wine2D, type = "p", pch = 19,
     col = df.prop$color,
     plot.type = "confplot",
     label.conf = c(TRUE, 1, 1),
     xlim = c(-0.9, 0.9),
     ylim = c(-0.9, 0.9))
abline(h = 0, col = "gray")
abline(v = 0, col = "gray")

par(pty = "m")

Subject Weights in 2D

Here we plot the weights of the individual subjects on each dimension. The code for extracting the weights from the results seems complicated but basically involves converting a list of matrices into a data frame based just on the diagonals of the matrices.

sswgts2D <- t(as.data.frame(lapply(wine2D$cweights, diag)))
par(pty = "s")
plot(D2 ~ D1, data = sswgts2D,
     xlim = c(0, 2),
     ylim = c(0, 2),
     main = "Subject weights in 2D")
abline(0, 1, col = "gray")

par(pty = "m")

Individual 2D Subject Spaces

Another plot we can make is the position of each stimulus in the individual subject space. These configurations are stored in the $conf attribute in the scaling results object.att (use the attributes(wine2D) command to see all the parts contained in the INDSCAL results object)

par(pty = "s")

par(pty = "m")

3D solution

wine3D <- smacofIndDiff(mlist,
                        ndim = 3,
                        type = "ordinal",
                        constraint = "indscal",
                        verbose = FALSE,
                        itmax = 5000)
wine3D
## 
## Call: smacofIndDiff(delta = mlist, ndim = 3, type = "ordinal", constraint = "indscal", 
##     verbose = FALSE, itmax = 5000)
## 
## Model: Three-way SMACOF 
## Number of objects: 10 
## Stress-1 value: 0.155876 
## Number of iterations: 1688

Plot the 3D solution without color

This plot cannot be printed since it is 3D and can be rotated by the user

plot3d(wine3D$gspace, type = "s",
       expand = 1.3,
       xlab = "X",
       ylab = "Y",
       zlab = "Z")
text3d(wine3D$gspace,
       text = df.prop$wine,
       adj = c(1.5, 1.5))

Plot the 3D solution with color

This plot cannot be printed since it is 3D and can be rotated by the user

plot3d(wine3D$gspace, type = "s",
       expand = 1.3,
       col = df.prop$color,
       xlab = "Wine Color",
       ylab = "Y",
       zlab = "Z")
text3d(wine3D$gspace,
       text = df.prop$wine,
       adj = c(1.5, 1.5))
sswgts2D <- t(as.data.frame(lapply(wine2D$cweights, diag)))
par(pty = "s")
plot(D2 ~ D1, data = sswgts2D,
     xlim = c(0, 2),
     ylim = c(0, 2))
abline(0, 1, col = "gray")

par(pty = "m")